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ABSTRACT 

We investigate whether tidal forcing can resuh in sound waves steepening into 
shocks at the surface of a star. To model the sound waves and shocks, wc consider 
adiabatic non-spherical perturbations of a Newtonian perfect fluid star. Because tidal 
forcing of sounds waves is naturally treated with linear theory, but the formation of 
shocks is necessarily nonlinear, we consider the perturbations in two regimes. In most 
of the interior, where tidal forcing dominates, we treat the perturbations as linear, 
while in a thin layer near the surface we treat them in full nonlinearity but in the 
approximation of plane symmetry, fixed gravitational field and a barotropic equation 
of state. Using a hodograph transformation, this nonlinear regime is also described 
by a linear equation. We show that the two regimes can be matched to give rise to a 
single mode equation which is linear but models nonlinearity in the outer layers. This 
can then be used to obtain an estimate for the critical mode amplitude at which a 
shock forms near the surface. As an application, we consider the tidal waves raised by 
the companion in an irrotational binary system in circular orbit. We find that shocks 
form at the same orbital separation where Roche lobe overflow occurs, and so shock 
formation is unlikely to occur. 



Key words: hydrodynamics 
- methods: analytical 



shock waves - stars:oscillations - stars: binaries: close 



X 



1 INTRODUCTION 

As far as we are aware it is unknown if the tidal forces in 
a binary inspiral can create shock waves before the binary 
objects touch, begin mass transfer or plunge. In order to 
investigate this, we have developed a quantitative criterion 
for the critical amplitude at which stellar perturbations form 
shocks that may be interesting in its own right, or in other 
applications. 

This work was originally motivated by the observation 
of unsmooth fluid behaviour in the numerical simulation of 
an irrotationa l, equa l mass neutron st a r (NS ) binary merger 
jBaiotti et all (|2008l 'l: iRezzoUa etlll (|2010l ')'). The simula- 
tions show surface waves breaking when the initial data are 
evolved with the cold equation of state (EOS) P = Kp^ , 
and a strong wind when they are evolved with the equivalent 
hot EOS P — pe (with initially constant entropy). As both 
simulations should be identical until genuine shocks form, it 
seems likely that both the wind and the surface waves are 
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artefacts of the interaction with the artiflcial atmosphere. 
On the other hand, these artefacts may also hide genuine 
shocks. 

Mass shedding in small amplitude nonlinear pertur- 
bations has been demonstrated numerically for neutron 
star s rotating near the ma ss-shedd ing angular v elocity 
(St crgioulas et all (|2006l '): [Pirnmclmc ver et al'l (|2004h j. but 
the minimum perturbation amplitude for this to occur has 
not been quantified. 

We have therefore tried to obtain a quantitative cri- 
terion for shock formation using a combination of stellar 
perturbation theory and nonlinear planar fluid dynamics. 
We consider a shock formation scenario where essentially 
radial sound waves steepen as they approach the surface 
because the density and sound speed approach zero at the 
surface. (Note that while the shape of the tidal bulges ro- 
tates around the star, individual fluid elements mainly move 
up and down.) 

If such waves are generated by tidal forces from the com- 
panion, their amplitude and shape is determined in the bulk 
of the star, where almost all the mass is. In this regime, lin- 
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ear perturbation theory can be used to obtain the response 
of the star to the tidal force, treating its proper oscillation 
modes as forced harmonic oscillators. For simplicity, we as- 
sume that the background star is irrotational and spherical. 

On the other hand, near the surface the fluid geometry 
can be approximated as plane-parallel, and entropy or com- 
position gradients become irrelevant compared to the den- 
sity gradient. In this regime we use a hodograph transform 
to cast the nonlinear dynamics into a single Imear second- 
order PDE. A shock forms if and only if the hodograph 
transform becomes singular: a criterion for this can be e x- 
amined within the model itself jGundlach fc Pleasj (|2009l )'). 
This criterion was tested in the numerical evolution of non- 
linear spherical perturbations of an n — 1 polytropic st ar 
and found to be accurate within 10% (|Gabler et al.l (|2009t )'). 
Once formed , these shocks quickly take a universal, self- 
similar form (|Gundlach fc LeVegu^ (|201ll )'). 

The two regimes are linked by noting that under certain 
approximations the fluid variables and their equations of 
motion in the two regimes coincide near the surface. 

In Sec. [5] we derive the combination of linear non- 
spherical and planar non-linear fluid motion. We summarise 
the (well-known) linear perturbation equations for adiabatic 
non-spherical stellar oscillations in a suitable notation, and 
their limit near the surface if the density vanishes there. We 
present the hodograph transform and the shock formation 
criterion. We cast the hodograph equation in a form where 
it can be identifled with the linear perturbation equations 
near the surface in a high-frequency Cowling approximation. 
We then evaluate the shock formation criterion on solutions 
of the linear perturbation equations as if they were solutions 
of the hodograph equation. 

In Sec. [3l we apply this general formalism to waves 
raised in a star by the tidal force of its binary companion. 
We can then use standard methods to calculate the reac- 
tion to this force by expanding the perturbations in proper 
oscillation modes. We obtain the orbital separation dcrit at 
which shocks first form as function of the modes of the star 
and the mass ratio q. 

In Sec. [4l we carry out the necessary numerical mode 
calculations for stars with polytropic equations of state. 
Sec. [S] reviews our main approximations and states our as- 
trophyical conclusions. 

A similar calculati on t o our S ec. [5] has been carried out 
for g-modes in NSs bv iLail (| 19941 ). Their frequency is lower 
than the orbital frequency at merger, and so the orbital fre- 
quency moves through resonance as the orbit shrinks, and 
the full time-dependent driven oscillator problem must be 
considered. It was assumed that no shock forms, and dissi- 
pative heating was estimated instead. It turns out that the 
duration of the resonance is too short to give rise to sig- 
nificant heating. By contrast, we focus on p- modes, which 
have higher frequencies and are never in resonance, and es- 
timate their amplitude adiabatically in the approximation 
of a stationary circular orbit. 



2 NONLINEAR EXTENSION OF LINEAR 
PERTURBATION MODES 

2.1 Linear adiabatic perturbation equations 

We consider linear adiabatic perturbations of a spherically 
symmetric static perfect fluid star in Newtonian physics, 
in the frequency domain. The background is assumed to 
be non-rotating and in hydrostatic equilibrium, d/dr is de- 
noted by a prime. The background quantities are the den- 
sity p(r), pressure P(r), gravitational potential 4>{r), grav- 
itational acceleration (p'{r) ^ 0, sound speed c(r) defined 
by = {dP/dp)s, entropy per rest mass s(r), and Brunt- 
Vaisala frequency N{r) defined by 



P_ 

P 



(1) 



The equations of hydrostatic equilibrium for the spherical 
background star are 



P' + 4>'p = 0, 
+ -0' = 47rG'/9. 



(2) 
(3) 



The displacement vector of the (polar) nonspherical 
perturbation is expanded in spherical harmonics as 

i{r,t) = e-"^' [Ur)Yim{e,ip)er + £.H{r)rVj_Yira{e,^)\ , (4) 

where 



V± = -ee 



d 



+ 



d 



T sin 6 dtp 



(5) 



Because the equations are linear, it is customary to make ^ 
complex as above for ease of calcuation. The physical dis- 
placement is its real part 



(6) 



We neglect axial displacements, which in a non-rotating star 
have no restoring force, and are zero modes. The (real) fluid 
velocity is simply the time derivative of the displacement, 
or 



oai(r,i) = Re[-ia;^(r,t)]. 



(7) 



The Lagrangian perturbation 5f of any background quantity 
/(r) is related to the Eulerian perturbation /i by 



5f=h+ i- V/ = /l + Cr/'. 



(8) 



All scalar perturbations are also expanded in spherical har- 
monics, for example the Eulerian density perturbation 

pi{r,t) = e''''pi{r)Yi^{e,^). (9) 

The assumption of adiabatic perturbations is that 



SP = cSp. 



(10) 



For given spherical harmonic index I, we define the short 
hand 

l{l + l)c^ 



F{r) 



1 



(11) 



The radial and (polar) horizontal parts of the Euler 
equation and the mass conservation equation can be com- 
bined to give a single second-order ODE for ^r- With a later 
approximation in mind, we write this as 



C" + AC + BC = S, 



(12) 
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where 

A = 

B = 



c F 



EL 

F 



N2_ 

9 



F 



<t>'i 1(1 + 1) 



T 



J F' 2 
-2- + — + - 
c b r 



(13) 



(14) 



(15) 



This is complemented by tlie perturbed Poisson equation 

(16) 



*i + -<^i ^ — 5-^01 = 47rG/9i. 



Its source term is the Eulerian density perturbation 



Pi 



P_ 

F 



/(^+ 1) 



<^1 



(17) 



Finally, the complete perturbation can be reconstructed us- 
ing 



a(r) 



Fu'^r 



+ oil 



(18) 



Our second-order equations for and 01 can be de- 
rived from t he fir st-orde r systems given in, for example, 
lUnno et alj (|l989l ) and IChristensen-Dalsgaard fc MullanI 
l| 19941) . 



2.2 Expansion near the surface and boundary 
conditions 



The boundary conditions for (jii are (jii 



at r = and 



0'i + + l)0i/'' = at r = R. The boundar y condition for 
at r- = is ~ r'"^ (|Unno et alJ (|l989l )'). To find the 
boundary condition for at r = _R, we need to expand the 
equations to leading order in x = r — i? 0. In the following, 
0{x) will be shorthand for 0(\x\/ R). 

We assume that near the surface the EOS is approxi- 
mated by the Gamma-law EOS 
pe 
n 



(19) 



where P is the pressure, p is the mass density, e the internal 
energy per rest mass, and n a constant. From the first law 
of thermodynamics, this is equivalent to 



(20) 



where s is the entropy per rest mass. The form of K{s) does 
not matter for our purposes (additional input is required 
to fix it), and more generally K can be considered a func- 
tion of both entropy and composition to take into account 
stratification effects in NSs. 

For a stratified background star where K{s) is a given 
function K{r), we find 



K' 



<P' n + 1 K' ^^^^ 

In the literature (for example IShapiro fc Teukolskvl (|2004h '). 
a power-law stratification is often considered which is pa- 
rameterised by specifying the equilibrium pressure as 



P(r) = Kop{r 



K{r) = Kop'^" ' 



(22) 



with no > n required for stability. (Note that n characterises 
the EOS, n/no characterises the stratification, and no ap- 
pears in the Lane-Emden equation). With p ~ x"" as a:: — > 0, 
we have 



Kir) 



(23) 



and so K, K' and N all diverge as x ~> 0. Clearly this 
assumption is unphysical on sufficiently small scales, and so 
near the surface the stratification must be adjusted away 
from a pure power law. 

To make an alternative quantitative assumption, we 
merely assume that A'''^ is finite at the surface, or 

iV^ = 0(l). (24) 

In this approximation, the condition of hydrostatic equilib- 
rium gives 

2 



c 

where 
9 = ^'(R) 



9^ I ,0/ 2^ 

\-0(x , 

n 



(25) 



(26) 



is the gravitational acceleration at the surface. It is also clear 
that 



0i=O(l), 0'i = O(l). 
We shall also need 
cb" = ^-g + 0{x), 



(27) 



(28) 



which follows directly from the background Poisson equation 
assuming n ^ 1. In this limit we have 

n + 1 



A = 
B = 



X 

n 

L 

gx 



+ o(i) 

1 



2(n-fl) 1(1 



1) 



+ o(i) 



gx 



n 
gx 
n 

gx' 



+ o(i) 



/(; + i) MR) 

na^ R 



+ 0(l) 



where 



.R 



a = LJ — — LJ 
9 



2_R1 

GM' 



(29) 



(30) 



(31) 



(32) 



is a dimensionless mode frequency. Note that uj and S axe 
defined as constants. 

Keeping only the leading 0{x~^) in A, B and S, (|12|l 
near the surface becomes 



d _^ n_ 

dx'^ : 



I d^r -2, n - 

; LO E,r = b. 

ax gx gx 



The solution of H33|) that is regular at a; = is 

ir{x) = C,U 



Anuj'-^x 



9 



S 

+ —, 

1 . 1-^ 



where we have defined the function 
/„(2) = 2"r(n-f 1) z-^J^{z). 



(33) 



(34) 



(35) 
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(This is a regular even function of z, normalised to obey 
/ii(0) = 1.) In the full set of equations S is of course not 
given a priori but is itself proportional to via (jlSp and 
H16p. We also need to fix an overall factor in the mode, and 
we choose to make the mode ^ dimensionless and set 



(36) 



Then C has a definite value (for any given mode and poly- 
tropic index n) , which is determined by solving the full equa- 
tions (|12[) and (|16|) . (In the Cowling approximation, where 
= 0, we would have C — 1.) 
The regular solution (|34|) obeys the boundary condition 



{n + l)g 



0. 



(37) 



T his b o undar y condition is equivalent to Eq. (17.69) 
of Cox (1 19801) ■ and the boundary co nditions derived in 
IChristensen-Dalsgaard fc MuUaril (|l994 ) but is not equiva- 
le nt to 5P/P = 0. Th is latter boundary condition is derived 
m lUnno et all (| 19891 ) under the assumption of finite sound 
speed at the surface, see their Eq. (18.31), and is therefore 
not applicable here. 

We introduce the dimensionless radius and mode fre- 



quency 
r 



(38) 



(Later, when we consider binaries, R and M will refer to 
Ri and Mi.) We can then write the approximation near the 
surface as 



Us) = Cfn (V4n<J'(s-l)) 



+ 



S 



(39) 



2.3 Nonlinear isentropic perturbations in the 
constant gravitational field, plane-parallel 
approximation 

In ICundlach fc Please! (|2009l ). we considered nonlinear 
smooth adiabatic motions in the approximations of planar 
geometry, the barotropic EOS 



(40) 



with K constant, and constant gravitational acceleration g, 
and derived the linear partial differential equation 



2n + l 
1 



in the independent variables 
fj. = 2nc, \ = V + gt. 



(41) 



(42) 



Here suffices denote partial derivatives. The surface is now 
at ^ = 0, and the interior of the star at > 0. The boundary 
condition = at = selects the regular solution. 

The criterion for the nonlinear fluid equations to form a 
shock is that the transformation from (x, t) to (n, X) becomes 
singular. We showed that this is equivalent to 



(1- 



- vx 



)' - v; < 0. 



(43) 



(Obviously, v must be real in this formula.) If and only if this 
condition is obeyed, a shock has formed, and the solution of 
pif) no longer has physical significance. 



2.4 Matching the two approximations 

We now have two sets of approximation: in the "pertur- 
bation approximation", everything is linearised around a 
spherical equilibrium solution. In the "hodograph approx- 
imation", the vertical fiuid motion is treated in full nonlin- 
earity, but we neglect horizontal motion, entropy gradients, 
and angular dependence, and approximate the gravitational 
field as fixed and constant in space and time. 

We expect that there is an overlap region just below the 
surface where both sets of approximations hold at the same 
time. In that region, we should then find the same equation 
of motion. To see this, note that the perturbation equations 
can be adapted to plane-parallel motion by formally setting 
I = and 1/R — (and hence F = 1 and a = 0), and 
to a constant gravitational field by setting ^'(r) = g and 
01 — 0. Neglecting entropy gradients corresponds to setting 
iV^ = and = —ngx (with g constant). With all these 
approximations, (|12[) reduces to 



d (,r n + 1 d^r 
rfx^ X dx 



LO^r ~ 0. 

gx 



(44) 



We now work from the other side. Consider a real, A- 
periodic solution of (|4ip of the form 



ii(A, /i) = Re (— ia;)e '"'^(x), 
where we have defined 



Ang' 



t ; 



A 



Then £_{x) obeys 



H —-fi rt^ C = 0, 

x'' X dx gx 



(45) 



(46) 



(47) 



which is of course formally the same equation as (|44p . al- 
though it represents nonlinear physics. Consider now a so- 
lution of (|4ip that represents a small perturbation about 
the hydrostatic equilibrium solution (|25p with n = 0, in the 
sense that 

\v\ < c, \5c] < c. (48) 

Then from the definitions (|42p and (|46p we can infer that 

x~x, i~t. (49) 

Furthermore, identifying the planar velocity v with the ra- 
dial velocity Vr , comparing (|45p with ([TJ , and using (|49p , we 
have 



(50) 



We have now justified the coincidence of (|44p and (|47p . 

However, the actual limit of the perturbation equations 
near the surface is not (|44p but (|33p . They differ in that to 
is not u, and by the (constant in x) source term S which is 
not present at all in the hodograph approximation. Tracing 
the differences back to the Euler and Poisson equations, we 
see that the terms in (|30l3ip proportional to l{l -\- 1) arise 
from horizontal motion, and the middle term in (|30p arises 
from the spherical (rather than planar) symmetry of the 
background. The difference between Cj and cj vanishes in 
the high-frequency limit cr^ S> 1. 

Generally, the source term S in (|12p represents the ef- 
fect of the perturbed gravitational potential on the fluid 
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displacement. In the hodograph approximation, such a term 
cannot be accounted for because the matliematics require 
the gravitational field g to be constant. However, the part 
S/ui'^ of the near-surface approximation (|39|) to the mode 
is constant in space, and so corresponds to the whole near- 
surface region bobbing up and down as one. Clearly, this 
part of motion has no effect on shock formation. We will 
therefore identify ^ with minus its constant-in-a; part. 

(We note in passing that in the high-frequency ap- 
proximation S/lj'^ ~ 4>'i{R)/^'^- If the mode oscillates with 
its own proper frequency, the corresponding displacement 
(j)'-^{R) /lo'^ cosujt is precisely what results from the gravita- 
tional field —(j)'i{R) cosojf.) 

In summary, to piece together the two approximations 
into a single "nonlinear mode equation" , we solve the stan- 
dard linear perturbation equations for and cpi on the 
whole domain ^ r ^ R, and then consider only the Cfn 
part of when we evaluate the shock formation criterion. 



Introducing the new shorthands 

^2 

^{n,K,y) = ^Uiyy + i^fUyr, (60) 

<I'(n, k) = maxtp{n, K,y), (61) 
y 

we can write out the shock formation criterion for a single 
mode as 

AcCacr^>I'(n,K<,)max|yi„| > 1, (62) 

Analysis of the function /„ shows that for k < 2{n + 1), 
i^iy) has its global maximum at y = 0, so ^' = 1. For k > 
2{n + l), the global maximum is attained at the point where 
fn + K/n = 0, with value > 1. In particular, if the mode 
oscillates at its proper frequency, ujf — cja, we have k = 1 
and hence Vl' = 1. 



2.5 Evaluating the shock formation criterion for a 
single mode with periodic time dependence 

Consider now a mode with proper frequency uja that 
is driven with another frequency cuf, resulting in some di- 
mensionless amplitude Aa. (To consider a mode oscillating 
freely, we can just set cj/ = oj^ in what follows.) Hence the 
actual time-dependent displacement radial displacement is 
given by 



Er{r,9,(p,t) ^ AaR C0S{UJ ft) Ylm{0,if) £,a{r-). 

Near the surface this approximates as 

Er{r,9,lfi,t) ~ AaRcOs(ujft) Ylr„(9,lfi) 



Col fn 



(51) 



(52) 



From the identification we have discussed, going into the 
bobbing frame we obtain 



V(X, /i) = -AaCaRujfYljn siu I —A ) /„ ( — ^ 



(53) 



Here we have, somewhat arbitrarily, chosen to use uja as an 
approximation to Ua, and (jj'iiR) an an approximation to 
S. In this formula, we consider w as a slowly varying function 
of the angles. 

Substituting (|53|l into the shock formation criterion 
and first focussing on the A-dependence, we can write 
the result as 



(1 - Ucosrf - (Vsinr)^ > 0, 
where 



y = — fJ; 

g 



9 



U = -A^C^Y,^ajf(y), 

V = AaCaYlm<^faaf'{y). 



(54) 



(55) 

(56) 
(57) 



We find that this criterion is sharpest periodically in r when 
COST = {//({/^ -I- V"^), and hence the criterion over at least 
a full oscillation period is equivalent to 



2.6 Evaluating the shock formation criterion for a 
single mode with arbitrary time dependence 

Consider now a mode driven with an arbitrary time- 
dependent amplitude "Zait), that is 

■Er{r,e,if,t) = R-E^{t) Yi^{e,ip) e„(r). (63) 
Hence, near the surface and in the bobbing frame 



V = RCa'B.a[^\ Yl„i{6,tfi) f„ (^fl 



(64) 



Hence we can write in dimensionless form 

vx ^p'{z)q{y), Vf,=p{z)q'{y), (65) 

where 

p{z) = (jJa^Ea (— ) , (66) 



My). 



(67) 



We then have to minimise (1 — vx)^ — vf^ over both z and y, 
to see if it reaches a negative value. 

Note that the periodic case of the previous subsection 
is recovered with p{z) = sinr, with r = ^/kz. 



3 PERTURBATIONS RAISED BY TIDAL 
FORCES 

3.1 Calculation of the tidal acceleration 

Consider a binary system with masses M\ and A/2 in an 
elliptic orbit. Let the orbital angular velocity and the spins 
of the stars with respect to an inertial reference system be 
n, 5*1 and 52- In the (non-inertial) reference system that 
moves and spins with star 1 and with origin in its centre of 
mass, the Euler equation for star 1 becomes 



+ [v- V)v + + V<?!-i = a, 

at p 

where 



(68) 



+ > 1. 



(58) 



a = —V(j)2 + ro — 25*1 X V — Si x f— Si x Si X f. 



(69) 
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where v, p and P are the fluid velocity, density and pressure 
in star 1, i^i and 4>2 are the gravitational potentials gener- 
ated by star 1 and star 2, respectively, fo{t) is the location 
of the centre of mass of the binary, and eo(t) the unit vector 
in its direction. Using angular momentum conservation in 
the form d\Q. = const, we find 



(70) 



In the following we assume Si = 0, both because this 
is beli e ved to be correct for NS binaries jBildsten fc Cutierl 
1 19921 ) : [Kochanei3 (|l992l )) and because it simplifies the cal- 
culation, as we can use perturbation theory on a spherical 
background star 1. The distance between the centre of mass 
of star 1 and the centre of mass of the binary is |ro(t)| = di 
and the distance between the centres of mass of the two stars 
is |r2(t)| = d = di + d2. Here Midi = M2d2 by virtue of fo 
being the centre of mass. 

We approximate 4>2 as spherically symmetric, that is 



and expand V(f>2 up to 0{r^) in f. Then 



(71) 



di + 



GAh 

d2 



diil 6*0 



GM2^_^,GM2^ ^ „, 2. 
-r + 3 — ;:^J-||o + 0{r ) 



d3 



d3 



(72) 



where _L Q, denotes the projection into the plane normal to 
Q and ||0 denotes the projection into the direction eb. The 
0(r*') term in round brackets vanishes by the assumption 
that the origin of r is the centre of mass of star 1. The 
remainder can be written as 



a = Vx + 0(r 



with 



GM2 2 SGM2 
' 2d^ ^ 2d3 



X COS ((> + y sin < 



(73) 



(74) 



where we have chosen Cartesian coordinates so that the orbit 
is in the xy-plane, and where (j){t) = J Qdt is the orbital 
phase, n{t) the instantaneous orbital angular velocity, and 
d{t) the instantaneous orbital separation. We can write the 
tidal potential in terms of spherical harmonics as 



GM2 2 , 

r + 

2d3 ^ 



^^^r^ sin^ 9[1 + cos 2(0 — ip)] 
4a^ 



3GM2 2 

r 

4d3 



sin^ 6 cos 2(0 — ip) — ^cos^ ^ \^ 



3GM2 2n 2^ , 
GM2 [W 2,, 



-2i0 



Y2, 



(75) 



Hence the tidal force is polar, and to leading order in r 
is given by quadrupole terms. To linear order in perturba- 
tion theory, the deformations caused by the Yqo, Y20 and 
y22 terms decouple. In circular orbits, the Y20 term is time- 
independent because d and Q. are constant, and therefore 
does not cause shocks. We can neglect it also for moderately 
eccentric orbits. The y22 term is always time-dependent be- 
cause (f> is: physically, the tides rotate around the star, so 
that individual fluid elements move up and down. 



3.2 Response of the star 

The response of the star to a is governed bv lLail (|l994h 

p^-f£^H(f,t)=pd, (76) 

where £ is a linear differential operator containing only spa- 
tial derivatives. 

It is generally assumed (for example, lUnno et al.l 

(|l989l )) that a star admits a complete set of eigenmodes 
which obey 



and 

where the inner product is 

(^'^ = Jf / pr^drdCl 



(77) 
(78) 

(79) 



1 

J [iimVi^ + I'il' + '^)^imViZ] pr^ dr. (80) 



In the second line, we have assumed a decomposition into 
spherical harmonics, and into the radial and horizontal 
parts defined by Eq. Q, and the standard normalisation 
J |y|^ dO, — 1. For simplicity of notation, and because a is 
polar, we neglect the axial parts of all vectors. Note that ^ 
is by our convention dimensionless, H has dimension length, 
and a has dimension acceleration. 

Decomposing (|76|) into modes, we have 



H(f,f) = i?^ Hc(t)el(r). 



where the dimensionless amplitudes 'Sa{t) obey 



_2 \ 1 {(,c,a) 



The general solution is 



Mt). 



(81) 



(82) 



Za{t) — BaS -I / smOJa{t — S) fa{s) ds. (83) 



3.3 Circular orbits 

For a circular orbit, only the I — m = 2 component of 
the tidal potential is time-dependent, through (fi{t) = 2Q,t, 
where Q is the constant orbital angular velocity. Neglecting 
the time-independent parts, and splitting the 3-dimensional 
gradient into a horizontal and radial part, we have 



GM2 



^3 . , ^Re e-^'"* (2ry22e- + r^V ^¥22 , 
d'^ y a 

Hence we have 

_ GM2 [6^ -2,nt (f.,(2s,s)) 
d^ V 5^ (e.,C.) ■ 

The solution of H76I) is then 



(84) 



(85) 



(86) 
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where the dimensionless amphtude of the particular integral 
is 



GAh 



67r {Cc,(2s,5)) 
5 



(87) 



Note that in contrast to lLail (| 19941 ). we approximate the or- 
bit as circular, and we do not take into account resonance. 
We assume that the tidal force is the dominant source of 
oscillations, and that the orbit evolves so slowly that tran- 
sients can be neglected, and so we set Ba = 0. 

We introduce the binary mass ratio and the dimension- 
less orbital separation 

Ma d 

We note 

G(l + q)Mr 



circ 

We can then write 



{Cc-,(2s,s)) 



A = J— g 

5 aW-Hl + q) 



(89) 



(90) 



Combining this with (|62|l , and using maxg \ Y22 \ ~ \J 15/327r 
and ijJf = 2fi, we have the dimensionless shock formation 
criterion 



q (C^,(2s,s)) 



> 1, (91) 



where a characterises the I 
the criterion, and where 



4(l + <?)' 



2 mode that maximises 



(92) 



3.4 Elliptic orbits 

For elliptic orbits, we have 



f22{t) = 



/20(t) 



fooit) 




l&K -2i.j.(t) (2S,S)) 

5 Cic.) ' 



d3{t) j V 5 

(Cc.(2s,s)) 
(el, el) 



^( ec,(2s,g)) 

{£,01, £,a) 



(93) 



(94) 



(95) 



For slightly elliptic orbits, the time-dependence of Q. and d 
is weak, and / = m = 2 is the dominant time-dependent 
deformation. For highly eccentric orbits, both and d~'^ 
are sharply peaked at periastron. This has two consequences: 
First, depending on how quickly perturbations set up by 
tidal forces are damped, it may be appropriate to treat each 
periastron passage as a transient, rather than as part of 
periodic excitation. Secondly, /22, /20 and /oo are now all 
comparably time-dependent. In fact, 



(96) 



where e is the eccentricity and dpa the periastron orbital 
separation. The three overlap integrals are also likely to be 
comparable. Moreover, as the orbital frequency even at pe- 
riastron is substantially lower than the mode frequency, we 



have approximately Ha {t) 



^ fa (t), as for the circular or- 



bit case. Hence we expect even the highly eccentric case to 
excite the 22, 20 and 00 perturbations at similar amplitudes 
'Ea{t), comparable to H22 in a circular orbit at periastron 
radius. 



3.5 Roche lobe overflow and resonance 

Our shock formation criterion becomes irrelevant once 
Roche lobe overflow occurs (or when two stars of identical 
mass and size touch). Defining 

dii 

where dn is the distance from the centre of star 1 to the 
Lagrange point Li, a simple calculation gives the quintic 

A"^ - 2A^ + A-'' - (1 + 3g)A^ + {2 + 3q)X - (1 + g) = 0. (98) 



(97) 



The one real solution gives Li (the four complex solutions 
give the other Lagrange points in the complex plane). Ap- 
proximating the Roche lobe as a sphere centered on star 1, 
Roche lobe overflow occurs at = A. As q S> 1, A ~ (3g)^^'^. 
Resonance is obtained at an orbital separation 



-2/3 (4^-4g)^/^ 



(99) 

Our result (|91|l has been obtained under the assumption that 
ujc, > 2Q, and this is always the case until shock formation 
or Roche lobe overflow. 



4 MODE RESULTS FOR POLYTROPIC STARS 

We have used a p ublicly available code 

(ICji ristensen-Dalsgaardl (|l997^ : [Christensen-Dalsgaard etiiil 
(| 19961 )) to calculate mode functions ^(s) and frequencies 
a and hence determine the overlap integrals and critical 
values of ?7 as a function of q. 

We have carried out the calculation for the isentropic 
(no = n) Gamma-law EOS 



P = Kp^ 



(100) 



with n = 1/2, 1, 3/2 and 3. n = 1 is often used as an 
approximate equations of state for cold NS matter, n = 
1/2 represents possible stiff neutron star equations 
of state, n = 3/2 and n = 3 are good approximations to 
non-relativistic and relativistic degenerate electron pressure, 
respectively. While these equations of state are simplistic, 
they have the advantage that the resulting stellar models, 
and hence our results, depend only on n, not on Mi and Ri. 
The mass, radiu s and polytropic constant are related by the 
scaling relation jShapiro fc Teukolskvl l|2004l )) 

(101) 



Ri cc K Ml 



In these simple stellar models there is no stratification and 
hence there are no g-modes. 

Our numerical results for circular orbits with equations 
of state n = 1/2, 1, 3/2 and 3 are shown in the figures. For 
the first several modes, we plot the critical value of orbital 
separation 77 against mass ratio. In all cases, shock forma- 
tion first occurs for the lowest frequency mode, so that is 
the curve that matters for shock formation. On the same 
plot, we also show the orbital separation -q at which Roche 
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0.1 1.0 10.0 100.0 0.1 1.0 10.0 100.0 

q q 



Figure 1. Critical value for shock formation of dimcnsionlcss 
orbital separation r) = d/Ri against mass ratio q = M2/M1. This 
applies for a circular orbit, the polytropic equation of state with 
n = 1/2 (r = 3), and I = m = 2 perturbations. The labels 
0, 1, 2, . . . denote the number of radial nodes in the mode. The 
critical value of r] for Roche lobe overflow through the Lagrange 
point Li is shown as a dashed line. At large q, rj oc q^^^ for all 
curves. 



10 




0.1 1.0 10.0 100.0 

q 



Figure 3. The same plots for the polytropic equation of state 
with n = 3/2 (F = 5/3). 



10 




0.1 1.0 10.0 100.0 

q 

Figure 4. The same plots for the polytropic equation of state 
with n = 3 (r = 4/3). (The bottom right corner appears blank 
because we have not calculated modes there. All our figures show 
the same ranges of rj and q.) 



Figure 2. The same plots for the polytropic equation of state 
with n = 1 (r = 2). 

lobe overflow starts. For all three equations of state, for all 
values of q, the critical orbital separation for the formation 
of tidal shocks coincides, within our approximations, with 
the critical orbital separation for Roche lobe overflow. 



5 CONCLUSIONS 

This paper consists of two parts: a quantitative criterion 
for linear perturbation modes to form shocks, and an ap- 
plication of this criterion to modes raised by tidal forces in 
compact binaries. 

Our shock formation criterion relies on the hodograph 
transformation to link linear perturbation modes in the in- 
terior to the fully nonlin ear shock formation criterion of 
iGundlach fc Please! (|2009l ) near the surface. This criterion 
is exact for plane-symmetric motion of a polytropic fluid in 
a constant gravitational field. The approximation of planar 
symmetry near the surface is natural, and it turns out that 



any buoyancy (non-barotropic) effects can also be safely ne- 
glected near the surface as long as the entropy gradient and 
any composition gradients are merely bounded. 

Our calculation of the tidal waves in perturbatio n 
theory is straightforward for circular orbits (|Lail (11994 )). 
For stars with a simple polytropic equation of state P = 
Kp^'^^^'^ in irrotational circular binary orbit, we find that 
the critical orbital separation for shock formation essentially 
coincides with the one for Roche lobe overfiow. In other 
words, tidal forces create shocks roughly when the 
binary begins to merge. Within our approximations, 
the two curves agree remarkably closely, so that we 
cannot say which actually occurs first. In any case, 
the p-mode shock formation mechanism we have in- 
vestigated here is not the primary mechanism for 
binary disruption. As discussed above in Sec. 13.41 
we expect the same result even for highly elliptic orbits. Al- 
though this is a negative result, it should be stressed that 
it was not obvious from dimensional analysis: we have also 
estimated the dimensionless factors. 

Extending our analysis to more realistic stellar 
models would require more extensive modelling, in 
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particular a realistic treatment of the surface. (We 
have shown above in Sec. 12.21 that another simple 
stellar model, assuming one polytropic constant for 
the equation of state and another for the stellar 
structure, gives rise to a divergent Brunt- Vaisala 
frequency, and so is inconsistent with our assump- 
tions of a perfect fluid surface). However, the fact 
that our result holds for a polytropic index ranging 
from n = 1/2 to n — 3 suggests that other equations of 
state would not show shock formation before merger 
either. Intuitively, the weakness of the shock forma- 
tion mechanism is dominated by a factor of (tidal 
force frequency /mode frequency)^ [the factor in 
Eq. (f62l)]. 

In a related result, iRosswog etal] (|2009l ) eive a cri- 
terion for the tidal disruption of a WD in orbit around a 
much more massive compact object 2 (black hole or NS) as 
(in our notation) 77disrupt — q^^^ based on numerical simu- 
lations. Our results are consistent with this for large q in 
that both Roche lobe overflow and tidal shock for- 
mation occur at the same 77, namely ryciit — lAAq^^^ 
for g » 1. 

Finally, this paper is motivated by the obser- 
vance of surface shocks in numerical simulations of 
binary neutron star merg er s .just before the sta rs 
touch (|Baiotti et al.l (|20oi ); iRezzolla et al.l (|2010h '). 
Our results indicate that these surface shocks are 
probably not physical. 
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